function S = salinityFunc_dev(P, T, C)
  % Function for calculating salinity ala SBE Data Procesing
  %
  % Author: David Goldsmith, Wash. State Dept. of Ecology, dgol461@ecy.wa.gov
  % Release date: N/A

  % Constants from IEEE Journal of Oceanic Engineering, Vol. OE-5, No. 1, January 1980, page 14
  % (via SeaBird Technical Note 14) 
  % Implementation verified 4/12/2012
    k = 0.0162;
    C35150 = 42.914;
    A = [2.070e-5 -6.370e-10 3.989e-15];
    B = [3.426e-2 4.464e-4 4.215e-1 -3.107e-3];
    G = [6.766097e-1 2.00564e-2 1.104259e-4 -6.9698e-7 1.0031e-9];
    a = [0.0080 -0.1692 25.3851 14.0941 -7.0261 2.7081];
    b = [0.0005 -0.0056 -0.0066 -0.0375 0.0636 -0.0144];
    S = zeros(size(T)); RP = S; RT = S;
  % Begin calculation
    R = C / C35150;
    val = 1 + B(1)*T + B(2)*T.*T + B(3)*R + B(4)*R.* T;
    good = (val~=0);
    RP(good) = 1 + (P(good) .* (A(1) + P(good) .* (A(2) + A(3)*P(good)))) ./ val(good);
    val = RP .* (G(1) + (T .* (G(2) + T .* (G(3) + T .* (G(4) + G(5)*T)))));
    good = (val~=0);
    RT(good) = R(good) ./ val(good);
    RT(RT <= 0) = 1e-6;
    tmp = RT';
    tmp = [ones(size(tmp)); tmp.^0.5; tmp; tmp.^1.5; tmp.^2; tmp.^2.5];
    sum1 = a * tmp;
    sum2 = b * tmp;
    val = 1.0 + k * (T - 15.0);
    good = (val~=0);
    S(good) = sum1(good)' + sum2(good)'.*(T(good) - 15.0)./val(good);
    S(~good) = NaN;
end